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Abstract 

We analyse the classical and quantum behaviour of a particle trapped 
in a diamond shaped billiard. We defined this billiard as a half sta- 
dium connected with a triangular billiard. A parameter £ which grad- 
ually change the shape of the billiard from a regular equilateral triangle 
(£ = 1) to a diamond (£ = 0) was used to control the transition be- 
tween the regular and chaotic regimes. The classical behaviour is regular 
when the control parameter £ is one; in contrast, the system is chaotic 
when £ / 1 even for values of £ close to one. The entropy grows fast 
as £ is decreased from 1 and the Lyapunov exponent remains positive 
for £ < 1. The Finite Difference Method was implemented in order to 
solve the quantum problem. The energy spectrum and eigenstates were 
numerically computed for different values of the control parameter. The 
nearest- neighbour spacing distribution is analysed as a function of £, find- 
ing a Poisson and a Gaussian Orthogonal Ensemble (GOE) distribution 
for regular and chaotic regimes respectively. Several scars and bouncing 
ball states are shown with their corresponding classical periodic orbits. 
Along the document the classical chaos identifiers are computed to show 
that system is chaotic. On the other hand, the quantum counterpart is 
in agreement with the Bohigas-Giannoni-Schmit conjecture and exhibits 
the standard features for chaotic billiard such as the scarring of the wave- 
function. 
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1 Introduction 

A billiard is a system where a particle, with mass m, is trapped in a region Q 
with perfect reflecting walls. The dynamics of the particle varies depending on 
the shape of the billiard boundary dD. The cardioid billiard, the Bunimovich 
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billiard (stadium billiard) and non-equilateral triangular billiards are typical 
examples which exhibit classical chaos. The quantum problem is reduced to 
solve the Hemholtz equation for the wave function ip(f) 

(V 2 + ft 2 ) ip(r) = for r G D (1) 

with the Dirichlet boundary condition xp(r) = if f G d D where ft = \j2mE /K 
is the wave vector and E is the energy. The energy level statistical properties of 
several Hamiltonian systems may be studied borrowing results from the random 
matrices theory. For example, it is well known that the energy level spacing dis- 
tributions of a system is Poissonian if its classical counterpart exhibits a regular 
motion. This is the case of billiards whose shape is a rectangle (particle in a two 
dimensional box), an equilateral triangle, a circle or an ellipse. On the other 
hand, if the classical counterpart has a chaotic motion, then the energy levels 
follow a Gaussian Orthogonal Ensemble (GOE) distribution pQ. Other non con- 
vex and chaotic two dimensional quantum cavities are the Sinai and the annular 
billiards. They introduce an inner disk of infinite potential into the rectangular 
and circular billiards respectively. 




Figure 1: (left) The diamond shaped billiard is defined through the functions 
f(x) and g(x). In turn, these functions depend on the parameters R, d\ and d^. 
(right) Triangular billiard. 



The system studied in this paper is the diamond shaped billiard (see Figure 
0. The upper boundary of the billiard is a half stadium defined by the equation 
V = f(x) with 

f{x)={ R ifM<f (2) 

i? 2 -(x-f) 2 iff <x <iL + R 
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and the lower boundary as y = g(x) with 

-d 2 x 



R 



di 

d 2 x 



-d 2 
-d 2 



if -R - \ < x < 



if < x < R + 



(3) 



The parameters which define the billiard depend on the control parameter £ 

according to: R(£) = ^(1-0, <M0 = (f+f)^ and d 2 (0 = ^fdi(0- 
The value of R Q has been taken as one and the non-dimensional parameter £ 
varies in the interval < £ < 1. The shape changes from a diamond to an 
equilateral triangle as the parameter £ goes from to 1. The boundary may be 
conveniently expressed in polar coordinates as follows 



r c (4>) 



R 

sin cf) 

r_(0) 



if </>a < ■ 
if < 
if 0c < 1 
if <\>d < 
if 6e < > 



< <I>b 

< 4>c 

» < 

< 2tt 



(4) 



where the left (— ) and right (+) quarter of circles are defined by 

r± (0) = ^ ^±di cos cj) + ^AR 2 -d\ sin 2 ^ . 
The lower points are located at 

d 2 



=pd2 cos 



— sin ( 



arctan 



and the ^-coordinate of the points from A to E are: (\>a = 0, 0, 
0c = arctan (^7) + 2 arctan (^), <j)D = ^ and = |tt, respectively. 



(5) 

(6) 
(*)■ 



2 Classical diamond billiard 

2.1 Trajectories 

There are two degrees of freedom in the diamond billiard, thus its phase space 
has four dimensions. Because of the conservation of energy the number of di- 
mensions is reduced to three. Commonly, in order to obtain the dynamical 
information of the system, we construct the Poincare section, and so we may 
study a two-dimensional map. This is equivalent to choose two variables which 
define where and how the collisions occur in the billiard. 

We chose a rescaled arclength l(<fr), and the angle a which defines the direction 
after the impact as the variables to describe the particle motion in the billiard. 
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Figure 2: Variables of the reduced phase space. The blue dashed line represents 
l((j>) while a is the angle formed by the normal vector at position </>, n(</>), and 
the trajectory of the particle after the impact. 



These variables are shown in Figure |2j The rescaled arclength is defined as 
l{4>) = j^k where L(0) = r c ((j))d(j). In order to compute the position and ve- 
locity of the (n+l)-th collision, using the position [x^ n \ y^) = (r c ((/>( n )) , </>( n )) 



and the incident velocity 



(vx\vy l) ) of the n-th collision, we may pro- 



ceed as follows: (i) Incident velocity of the (n+1) collision. The normal and 
tangent v[ n ^ components of the velocity are computed by projecting it into 
the normal, n((/>), and the tangent, £(0), unitary vectors of the boundary. Just 
after the collision with the boundary the normal component velocity changes 
its sign while the tangent component remains unchanged, thus one can obtain 
the velocity ^ n+1 ) after the n-th collision which is also the incident velocity 
(n + l)-th collision. For this calculation the components of the tangent vector 
t = jt (r c ((j)) cos (/>, r c ((j)) s'mcj)) = (T x: T y ) are needed. These are 



and 



T x (4>) = < 



— r+(0) sin cj) + cos (j)w+((j)) if (j>A < « 

R if (j) B < 

r -((/)) sin^ + cos (j)W-((j)) if <pc < 1 

{R + d 1 /2)a if (j) D < 



where 



and 



r+(0) cos^ + sin (j)w + ((j)) if <pA < <j> 

if <j> B < <j> 

r y ((^) = \ t -((/)) coscj) + sin (j)W-((j)) if <pc < (j) 

—d2d if <j)]j < <j) 

d^a if <Pe < 4> 



(A\ ^ l J • A fl -L ^1COS0 

™±(0) = T^^i sm0 1 ± - — - 

2 \ 2r±(0) =p di cos0 



< 

' < 4>c 

< <\>D 

> < 2tt 
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< <\>D 

< (t>E 

< 27T 



(R- 
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(8) 

(9) 
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The normal vector is obtained by rotating the tangent vector n 



-tt/2)£, 
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hence n x ((j)) = t y ((j)) and n y ((j)) = —t x {<j>)- (w) Position of the (n+l)-th collision. 
If the line which crosses through the points {x^ n \y^ n ^ and (x^ n+1 \ 2/( n+1 )) is 

m {n) x . 



Y^ n \x) 



, then the slope and the ^/-intercept are 



rn 



(n) 



(n) 



(n) 



in) 



. , and b M =y 

(n) * (n) e 



An) 
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respectively. The intersections of a line Y(x) = mx + b with the boundary are 

R and Y(x%) > 



m 
X_ 

X + 

X- 



if f < x% < f 
if 1^1 < f andF(^) >0 



di 



< £± < and F(x^) > 



if — R 2 ^ du — — 2 
if -i? - ^ < x + < and F(x + ) < 
if < x_ < R + ^ and < 



(12) 



where we have defined 



(b + d 2 )(di + 2#) 
■2i?) 



2c?2 i m(di 
± _ s di - 26m ± ^4(1 - 



and 



2 )i? 2 - (dim + 26) 2 



2(1 + ra 2 ) 



(13) 
(14) 



with s = {+, — }. The diamond billiard is a convex billiard, so the equation (12) 
gives us two roots: x* and x\, one of them is the position of the current collision 
so it is known, let us call it x n = (xl)^ n \ the other root gives the position of 
the (n + 1) collision 



= (x*)^ and ?/ n+1 ) = ((x$) (n) ) . 



(15) 



In general this procedure works well. Nonetheless, if the particle reaches one of 
the points {A, £?, C, D, E} where the tangent and normal vectors to the bound- 
ary are not defined, then the method fails. Although, this situation for an 
arbitrary initial condition rarely happens, the problem sometimes is solved by 
taking the average of the normal and tangent vectors for the boundaries con- 
nected in those problematic points. In Figure [3] are shown some trajectories for 
the triangular billiard and the diamond billiard using the procedure described 
above. 

2.1.1 Entropy 

In the previous section, a methodology based on geometry was used to find the 
classical trajectories of the particle in a diamond billiard. Indeed, this is not the 
more elegant way to find trajectories, and there should exist a transformation 
or map which connects the variables of the reduced phase space of consecutive 
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Figure 3: Several trajectories after 1000 collisions are plotted. (Upper panel) £ = 
1 — 1 x 10 -16 , note that regular behaviour is obtained only for values of £ very close 
to 1. (Middle panel) Some closed orbits of the triangular billiard. (Lower panel) For £ 
far enough from 1 (diamond billiards), the trajectory fills the billiard irregularly (left) 
f = 0.99 and (right) f = 0. 

collisions of the diamond billiard. In principle, the trajectories of the particle 
may be constructed with the knowledge of the billiard map and the initial 
conditions. If we avoid the very special cases of the periodic orbits, then the 
degree of irregularity of a set of trajectories with different initial conditions 
should depend only on the shape of the billiard. The entropy S is calculated 
in order to determine quantitatively such degree of irregularity. S may be 
computed as follows: Let a n be the incident angle with respect the normal vector 
on the boundary. The range of this variable is the interval I = [— 7r/2, 7r/2]. 
This interval is divided in M equal subintervals Ij. Then N collisions and their 
respective incident angles a n are generated, where n = 1, 2, 3, .., N. If Ni is the 
number of angles a n which live in the interval i^, then the probability to find 
an incident angle in the interval Ii is P (1%) = and the entropy S may be 
computed in the standard way as 

N 

S = -Y>(J;)ln[P (/;)]• (16) 
i=l 

The maximum entropy is obtained when the set of generated incident angles 
{a n } are uniformly distributed in the subintervals {Ii}. For this case the prob- 
ability is equal for each subinterval, hence P (Ii) = ^, and the entropy takes 
its maximum value S max = ln(iV). On the other hand, if all incident angles 
lie in a single subinterval 2j, then the probability would be P (Ii) = and 
the entropy is zero. The entropy computation using the equation ( p~6] ) generally 
depends on the initial condition used. In order to avoid this dependence we 
have computed the entropy for 1000 trajectories with different random initial 
conditions. Posteriorly, these entropies were averaged for each particular value 
of the control parameter (Figure [4]). The smallest value of entropy is obtained 
when £ is exactly one and the billiard is an equilateral triangle. The entropy 
grows quickly as the half of stadium is introduced in one of the triangle sides, 
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Figure 4: Entropy, (left) Entropy computed for a single value of £ and one initial 
random condition. The number of collisions of these graphs were N = 3000, and 
the interval subdivisions were M = 100. (right) We have generated 1000 random 
initial conditions with their respective entropies, and later these entropies were 
averaged. 

even when the control parameter is close to one as £ = 0.99 where the corre- 
sponding entropy is about the 55% of its maximum theoretical value. As £ is set 
far from one, the entropy practically stabilizes its value reaching about a 70% 
of S max . In this regime, the trajectory of the particle is more complex than the 
one found for £ close to one as it is clear from a comparison between the lower 
and upper panels of Figure [3] 

2.1.2 Lyapunov exponent 

The Lyapunov exponent A is used as a measure of divergence between trajecto- 
ries for a couple of infinitesimal close initial conditions in the phase space. The 
time is not a suitable parameter in order to compute A in billiard systems since 
the particle movement is linear while it does not collide and the trajectories will 
diverge linearly. The collision index n was used as parameter instead of time, 
as usual a great sensibility with small changes of the initial conditions is char- 
acterized by S n — S exp (An) where S n is the absolute value of the difference 
between incident angles of nearby trajectories after n collisions. 

The Lyapunov exponent was averaged in order to avoid the dependence with 
the initial conditions. For a couple of close trajectories random initial condi- 
tions were generated, and the Lyapunov exponent was computed for each initial 
condition (Figure |5j. The Lyapunov exponents computed in the previous step 
were averaged for each value of £ (Figure [6]). In order to minimize the error 
introduced by the saturation, we decided to calculate the slope between adja- 
cent points and average it for each single Lyapunov exponent computed, thus 
the saturated points frequently have small contribution due to the alternation 
of the slope sign. Some graphics are not well defined as the one shown in Figure 
[6] and some inaccuracies persist in the final result, even if the number of initial 
conditions is increased. For this reason, the final averaged on the Lyapunov 
exponent in Figure [6] is not as smooth as the entropy of the Figure [5] How- 
ever, the Figure [6] is able to capture an important feature of the billiard: as the 
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Figure 5: Lyapunov exponent saturation. This is a typical graphic obtained 
for a single value of the control parameter (£ = 0.9999) and one random initial 
condition; 40 collisions are considered. For this particular situation the satu- 
ration occurs about after 13 collisions. The Lyapunov exponent is the slope of 
the non-saturated part. 




0.0 0.2 OA 0.6 0.8 I.C Q.O 0.2 0.4 0.6 0.8 1.0 



Figure 6: Averaged Lyapunov exponent. The number of random initial condi- 
tions of each graph was (left) 10 and (right) 1000. 

half stadium appears over one side of the original equilateral triangle, then the 
Lyapunov exponent substantially increases, and the non negative values of it 
ensures a great sensibility to the initial conditions, even for values of £ close to 
one. 

3 Quantum diamond billiard: Finite Difference 
Method Implementation 

This problem is typically solved by using finite element method (FEM) and 
iterative methods on the discretized version of the Schrodinger equation [2] . We 
may use the FDM to express the Hamiltonian as a matrix on a lattice, and 
then solve the resulting eigenvalue problem. The discretization of the region is 
shown in Figure [7]- left, each point at the position is labelled in one of these 
ways: with pair or the single index u. The second option is used in order 
to avoid the impractical use of four indices in the Hamiltonian matrix. During 
the lattice construction it is easy to build the function u = u(i,j) which maps 
from the pair of indices to the point u, let us call this process as the first 
indexing. The time independent Schrodinger equation is evaluated at the point 
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Figure 7: First and second indexing, (left) The domain is discretized in a 
rectangular lattice of M x N points including the boundary (the points out the 
solid rectangle). Each point (i,j) is labelled by the single index u. (right) Here 
only the inner points are indexed, we build the map a = a(i,j), inversely the 
pair (i, j) is obtained by i = 12(0) and j = .72 (a)- The points of the two indexing 
are related by u = u(i,j) = u (^(a), h( a ))- The total number of inner points 
is Q = M' x N' 



u(i,j) according to 



2m 



v(m^)\ u(i , j) = Em\ u{i , j) - 



(17) 



The second derivatives of the laplacian may be evaluated using central differ- 
ences [3] 



d 2 iP 
dx 2 



and 



d 2 *P 
dy 2 



u(i,j) 



u(i,j) 



5x 2 



1 

Sy 



2 [^n(i,j + l) 



u(i,j-l) 



The notation is simplified defining 

=u(i(u) ± l,j(u)) . 



(18) 



This transformation makes a horizontal displacement in the lattice from the 
point u. Although, there are several values of u for a single i or j (TV values 
for the index i(u), and M values for the index j(u)) we have only one value 
for u fixing both i and j. The equation (18) gives the neighbor of u at its left 
(— ) or right (+). Similarly, the vertical displacement from the point u(i,j) is 
computed with 

u±=u(i(u),j(u)±l) . (19) 
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With this notation, equation (1171) may be written as 



^ Huyipy = E^ u where 



2m 



Sx 2 



Sy 2 



25 u , v (5x- 2 + Sy- 2 ) 



(20) 

is the Hamiltonian (repeated indices in the last term do not involve sum over 
them). Since the problem is solved only for the inner points, we performed a 
second indexing (see Figure fright), so the eigenvalue equation may be written 

as 

Q=M'xN' 
veD [3=1 

The eigenvalues and eigenvectors of the Hamiltonian H a p are obtained numer- 
ically. Commonly, the packages of matrix diagonalization arrange the eigenvec- 
tors in a matrix, let us call it M a p 



M = 



( Mn Mia 
M 2i M 22 



Miq \ 

M 2Q 



J 



(22) 



V Mqi Mq 2 • • • Mq 

If the eigenvectors are arranged in the columns of such matrix, then the state s 
evaluated at the point a is ipa^ = M as with (s = 1, 2, Q). In order to return 
to the initial labelling we may write 







if fij e D 
otherwise. 



(23) 



The generalization for billiards of arbitrary shape does not represent consider- 
able difficulties. We may place the boundary of the arbitrary shape billiard over 
the rectangular grid and take only the points inside of it. After the identification 
of the boundary points, the inner points (say Q inner points) are enumerated 
first (u = 1,2,...,Q), and the boundary points later, so a second indexation is 
avoided. Some grids for the billiard consider in this study are shown in Figure 

El 

The numerical results were compared with the exact ones for the equilateral 
triangular billiard. The triangular billiard is integrable and the expression for 
the energy levels is well known [4] 



E n — E n 



/4tt\ 2 / h 2 
V 3 J \2mdj 



(P 2 



pq) 



(24) 



where d\ is the edge length when £ = 1, p and q are positive integers which 
satisfy q G [l,p/2]. A comparison of the numerical and analytic energy levels 



10 



Figure 8: The grid. As the density of points is incremented, the mesh is better 
adjusted to the geometry of the billiard. Some grids are shown, the number of 
inner points are: (left) Q = 32, (center) Q = 1108 and (right) Q = 4618. 



Table 1: Triangular billiard eigenvalues. 



State 


E n (J) 


E n (J) 


Relative Error 


n 


FDM 


Exact 


% 


1 


2.14940 


2.14849 


0.04235 


2 


5.01655 


5.01313 


0.06817 


3 


5.01743 


5.01313 


0.08570 


4 


8.59995 


8.59394 


0.06988 


5 


9.31609 


9.31010 


0.06429 


6 


9.31662 


9.31010 


0.06998 


7 


13.6095 


13.6071 


0.01763 


8 


13.6171 


13.6071 


0.07343 


9 


15.0451 


15.0394 


0.03788 


10 


15.0451 


15.0394 


0.03788 


11 


19.3386 


19.3364 


0.01137 


12 


20.0559 


20.0525 


0.01695 


13 


20.0559 


20.0525 


0.01695 


14 


22.1997 


22.2010 


0.00585 


15 


22.2001 


22.2010 


0.00405 



was done and listed in Table 1. One advantage of the FDM lies in the fact that 
Hamiltonian is computed by a direct evaluation of the potential and some Kro- 
necker deltas, so in a personal computer (in our case an i7 processor) building a 
Hamiltonian matrix of 11000 x 11000 take less than a minute, and its orthogo- 
nalization with the lapack package using Fortran, requires about 25 minutes. In 
order to check the accuracy of the results, a comparison with the energy stair- 
case function Af(E) with the Weyl-type formula was performed. Af(E) gives 
the number of energy levels under the energy E and it is defined by 

M{E) :=J20(E-Ei) (25) 

i 

where 0(x) is the step function. The analytical result for a two-dimensional 
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billiard with area A and perimeter P is given by [5j [H El IE] 

Carefully speaking, the Weyl formula is valid in the semiclassical limit, that is for 
high energy levels. Other typical problems in billiards as the scar identification 



Figure 9: Diamond shape billiard wavefunctions. Here some eigenfunctions are 
plotted using the a grid with Q = 4618 inner points (see Figure [8]), and the 
control parameter was set as £ = 1. We have set m = 1 and H = 1. 

deep in the semiclassical limit may be faced using more convenient methods. 
For instance, the improved Heller's PWDM (Heller's plane wave decomposition 
method) [9]. This method avoids the computation of the eigenvectors near to 
the ground state and goes directly for the computation of the states with high 
quantum numbers. We used the finite difference method in order to diagonalize 
quantum diamond billiard. Some of the first excited states of this billiard are 
shown in Figure [9] In Fi gure [10| we have superposed the numerical result of 



the energy staircase function with the one provided by (26). The deviation of 
the numerical results to the theoretical prediction for high energies is due to 
numerical errors because the discretization procedure is not able to describe 
properly wavefunctions corresponding to very high energies, which have very 
small wavelength oscillations. From this, it is clear that only the first ~ 200 
computed states are reliable. 



4 Quantum diamond billiard: level statistics 

Several experiments were performed in the nineties with quantum hard wall 
billiards, e.g. the microwave resonators, which used the equivalence between 
the stationary Scrodinger equation with the Helmholtz equation to study chaos 
in quantum billiards using electromagnetic waves [10] . Other devices used in 
the quantum chaos study were the semiconductors billiards, those are open 
quantum cavities which permit a current flow through two contact points. These 
structures are different from a quantum billiard, which is completely closed and 
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Figure 10: Spectral staircase function for the diamond shaped billiard, (left) 
Superposition of the numerical staircase energy staircase function (solid line), 
with the Weyl's formula (dashed line), (right) A zoomed region corresponding 
to the blue rectangle of the graphic on the left side. 



confines a single particle inside it. However, if the size of the quantum open 
cavity is much smaller than the mean free path of the electrons, then the device 
approaches a quantum billiard and it shares with the quantum billiards several 
properties e.g. energy level statistics and the scarring of the wavefunction [2]. 
For Hamiltonian systems such as the one described in this writing the energy 
level spacing distribution, P(s), is a feature which distinguishes the spectrum of 
a system with regular or chaotic classical analogue. According to the Bohigas- 
Giannoni-Schmit conjecture [IT] the spectra of a time-reversal-invariant system 
with a classical chaotic counterpart follows a Gaussian Orthogonal Ensamble 
(GOE) distribution. On the other hand, if the classical analogue is regular, 
then the spectrum is characterized by a Poisson distribution (see Figure 11). 
This conjecture has been tested in a variety of systems including billiards. 




Figure 11: Nearest neighbor spacing distribution of the energy levels of the 
triangular billiard. The histogram was built using the analytic expression for 
the energy levels of the equilateral triangular billiard (see equation (24)). The 
level statistics is Poissonian (solid line) because the equilateral triangle billiard 
is regular. For the irregular case, the nearest neighbor spacing distribution will 
follow either a GOE2 distribution (dashed line) or a GOE distribution (dot 
dashed line) according to the billiard symmetries. 



The diamond shaped billiard has a mirror reflection symmetry axis. For this 
reason, there are two set of states related to each symmetry class, namely, odd 
or even eigenstates. The general expression of the nearest neighbour spacing 
distribution for a superposition of N independent spectra in the GOE statistics 
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is given by [T2] 



N GOE 



to 



r 

ds 2 



erfc 



~2N 



n AT 



(27) 



where s is the energy spacing between nearest neighbour levels, and erfc(x) is 
the complementary error function. The spacing distribution for N = 2 is 



^2 goe 0*) = 



7TS 



e 10 erfc 



/7TS 



(28) 



| 1 < state <800 g =0 | 




Figure 12: Nearest neighbor spacing distribution of the complete billiard, (left) 
The histogram of the level spacing distribution was built using the first 800 
energy levels, the first one hundred states were not taken. The total number 
of energy levels computed using the finite difference method was Q = 11026. 
(right) Ninth state of the complete diamond shaped billiard. 



In Figure 12 it is shown the nearest neighbor spacing distribution of the 



diamond shaped billiard which fits the P2GOE (GOE2) distribution, as expected. 

There are two ways to recover the GOE distribution: the first one is to 
classify the energy levels according to the parity of the eigenstates and the his- 
togram is built with one of the two sets. However, there is a disadvantage in this 
method, because it requires to take approximately the half of the energy levels 
computed. The second one consists in statistical study of the corresponding 
desymmetrized billiard spectrum. In this case the billiard is desymmetrized by 
taking a half of it for the mirror symmetry of the diamond billiard. The level 
statistic effectively obeys a GOE distribution for the desymmetrized billiard. 



The result is shown in Figure 13 



Another important feature of quantum chaotic billiards is the high concen- 
tration of the wavefunction amplitude along the classical periodic orbits. The 
phenomenon was initially observed by McDonald and Kaufman [13 j, and pos- 
teriorly in the Bunimovich billiard by Heller [14], who called it a scar. Using 
the analytic solution of the wavefunction it is not possible to built a scar in the 
rectangular, circular or equilateral triangle billiard. The scarring of wavefunc- 
tion in billiards does not appear in regular billiards, and it is exclusive for the 
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Figure 13: Nearest neighbor spacing distribution of the desymmetrized dia- 
mond shaped billiard, (left) The histogram of the level spacing distribution 
was built using the first 500 energy levels. The total of energy levels computed 
using the finite difference method was Q = 11252. (right) Ninth state of the 
desymmetrized diamond shaped billiard. 



chaotic ones. As the quantum numbers are increased, we expect to recover the 
classical characteristics of the system, which is, in some sense, the idea behind 
the correspondence principle. In a chaotic billiard, the trajectories which evolve 
from an arbitrary initial condition tend to fill the whole billiard, as consequence 
a typical wavefunction should not have a significant localization, which is the 
more common situation for the irregular billiards. Nonetheless, for the special 
case of an unstable periodic orbit, it is possible to find a high probability density 
underlying such classical trajectory, as we may intuitively expect at least in the 



semiclassical limit. Some scars and bouncing ball states are shown in Figure 14 
The bouncing ball states have a well defined momentum, but not position and 
we may associate a set of classical periodic orbits to a single bouncing ball state. 
In contrast, a scar is related to a single unstable periodic orbit. 



Example of scars are shown in Figure 14 -(a) and 14 "(c). The stability of the 



orbits with lowest period is shown in Figure [151 



5 Quantum diamond billiard: time evolution of 
the state vector 

In this section the study will be limited to the time evolution of the mean values. 
The time evolution of the state vector for quantum Hamiltonian systems is 
well known | — e~ % ^ | \I/(0)). Expanding the state vector over the 

corresponding stationary states | \£ (r, 0)) = ^ s c s \ ip^), and evaluating on 
one arbitrary inner point of the lattice we find 

(?\*(t))\r=r u = *«(*) = I> eX P (-^)^ } - 

The components are given in the usual way, however we may use the lattice in 
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Figure 15: Orbits stability, (left) Unstable periodic orbit after 35 collisions, 
(right) stable periodic orbit after 1000 collisions. 



order to compute them easily 



r,= lim y#$ B (0)^ (29) 



Since the eigenvectors of the Hamiltonian matrix, given by the equation (20), 
are real, then the complex conjugation has been dropped. If the number of 
inner points is large, we may use the last equation without the limit as a good 
approximation for the components computation. The same approach may be 
used for the computation of the several mean values involved in the uncertainty 
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products for momentum and position. Using the index displacement transfor- 
mations and central differences, the gradient evaluated at the point u(i,j) may 
be expressed as 



1 /¥„+(*)-¥„-(<) ¥ u+ (t)-¥„_(t) 



Sx 



Sy 



so the mean value for momentum takes the form 

h 



$(t)) = y lim Y, *«(*)* (^*) 



<5 2 f. 



(30) 



(31) 



Taking the real and imaginary part of the last equation we find 
$(t)) = h lim V [v,*l (t)S 2 r 



and 



lim V (V,*) (t)S 2 r = 



(32) 



where the following expressions were defined 
[V, *] (t) := Re [*„(*)] Im [(w) (*)] - Im [*„(*)] Re [(w) (t) 

{V, (t) := Re [*„(*)] Re [(w) (t)] + Im [*„(*)] Im [(v*) (t) 



and 



(33) 



The condition in the equation (32) appears because the mean value is a real 
quantity, then for practical means, the imaginary part may be used to check the 
accuracy of the numerical integral evaluation. Similarly, for the mean value of 
the squared momentum we find 



with 



(f 2 (t)) = -h 2 lim V{v 2 ,*} (t)S 

lie® 

lim \^ 2 ^} (t)S 2 r = 



(34) 



where the laplacian at the point u(i,j) is 



•2(fe- 2 + ^- 2 )^ n (t). (35) 



Finally, the average of an arbitrary function / with only position dependence is 
(m)(t) = lim V | *„(*) | 2 /„. (36) 

o 2 r^0 z — * 
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Figure 16: Gaussian wave packet evolution. (Top), (left) Initial state, (right) 
Uncertainty products: AxAp x (solid line) and AyAp y (dashed line). The dot 
dashed line is the minimum uncertainty value AxAp x = ft/2 (dot-dashed line). 
(Bottom). Trajectory of the position operator expectation values. 



The time evolution of the position mean values of the system prepared in a 
well localized initial state at f Q = (x ,y ) using a Gaussian wave packet 

- - (y^- +i( Ka x+ Ky y) 

(r | *(0)> = X ——^ (37) 

Z7Ta x a y 



is shown in Figure 16 k = (K x ,K y ) is the wave vector, a x and a y are the 
standard deviation on along x or y respectively. As seen in Figure [lT| the 
wavepacket is destroyed after a few collisions. However, this is not a conse- 
quence of the classical chaotic behaviour and such irregular dynamics may be 
attributed to a non-coherent preparation of the initial state. This is checked 
in the evolution of the uncertainty products. The system must be prepared in 
a coherent state in order to reduce the uncertainty products to their minimum 
value ft/2; nevertheless, we do not have a general analytical expression for the 
coherent states of billiards, even in simple cases such as a particle in a rectan- 
gular box. 

Classical chaos of an specific system often emerges from its nonlinear nature. 
Nevertheless, the classical and quantum billiard systems are an exception of this 
rule because of the absence of nonlinear terms in their governing equations. In- 
deed, the difficulty for quantum chaos determination does not lie in this lack of 
nonlinearity rather the problem resides in the difficulty to find a correspondence 
between the classical and quantum behaviour far from the classical limit when 
the classical system evolves chaotically. The question is not solved by simply 
proving the Bohigas-Giannoni-Schmit conjecture because the nearest neighbor 
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Figure 17: Probability density distribution of the billiard prepared in a Gaus- 
sian wave packet. The packet was initially placed in f Q = (— R, 0). The wave 
vector was pointed to 7r/4 rad respect the x-axis. Classically, the particle should 
describe the 4-periodic orbit shown in Figure [15] nonetheless, after the second 
collision the corresponding state is highly delocalized. The particle tends to 



remain on a bouncing ball state after the second collision (see Figure 14 -(f)); 
however, the next collisions increase the derealization and the evolution turns 
non-coherent. 



spacing distribution is just a semiclassical result. The analysis of the position 
operator expectation value is an alternative to study the correspondence be- 
tween the classical and quantum system in an irregular regime. However, this 
approach frequently is not successful because the quantum system evolves in a 
non-coherent way when its classical counterpart is chaotic as we show in this 
writing. This would be reason, for which the quantum Hamiltonian system sen- 
sibility has been sometimes studied by perturbing the Hamiltonian instead of 
by changing the initial state [15] . 



6 Concluding remarks 

The classical and quantum diamond billiard was studied through some quanti- 
ties. In particular, we calculate the entropy, the Lyapunov exponent and some 
trajectories for the classical problem. The classical chaotic behaviour emerges 
fast with small modifications on the boundary of the regular equilateral billiard 
(£ = 1). The entropy and the Lyapunov exponent grow when a half of stadium 
is introduced in one side of the triangular billiard. If the control parameter 
is set far enough from one, say in the interval 0.8 < £ < 1, then the entropy 
practically is a 70% of S max . This percentage is relative far from its maximum 
and it occurs because the diamond billiard does not have dispersive frontiers as 
other billiards e.g. Sinai billiard. Nonetheless, it is enough to ensure a great 
irregularity in the classical trajectories when the entropy is about 0.7S max . 
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The finite difference method provides a way to solve the quantum problem. 
The diamond shape billiard has a mirror reflecting symmetry. Because of this, 
the energy levels split in two different symmetry classes according to the wave- 
function parity. As consequence, P(s) for the complete billiard is given by a 
GOE2 distribution. If diamond billiard is desymmetrized, then the level statis- 
tics follows a GOE distribution. On the other hand, the classical behaviour is 
regular when the control parameter is set to one and the distribution is Pois- 
sonian. Therefore, the results are according to the Bohigas-Giannoni-Schmit 
conjecture. We found scarred states in the quantum diamond billiard, as well 
as bouncing ball states with their corresponding set of classical stable periodic 
orbits. These results are in agreement with previous work in the field for other 
Hamiltonian systems. In the last section, a practical way to compute the time 
evolution of the state vector was described and the lattice previously built in 
the finite difference method implementation was used for this aim. 
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